arXiv:1505.06142vl [math.DS] 22 May 2015 


On the N -Extended Enler System I. 
Generalized Jacobi Elliptic Functions 

S. Ferrer, F. Crespo, and F. J. Molcro 

Dpto de Matematica Aplicada, Universidad de Murcia, 30071 Espinardo, Spain 


May 25, 2015 


Abstract 

We study the integrable system of first order differential equations = 

a, Ilj^i u} j( v )i (If; hj < AT) as an initial value problem, with real coeffi¬ 
cients oti and initial conditions w,;(0). The analysis is based on its quadratic 
first integrals. For each dimension A, the system defines a family of func¬ 
tions, generically hyperelliptic functions. When A = 3, this system gener¬ 
alizes the classic Euler system for the reduced flow of the free rigid body 
problem, thus we call it A-extended Euler system (A-EES). In this Part I 
the cases A = 4 and A = 5 are studied, generalizing Jacobi elliptic functions 
which are defined as a 3-EES. Taking into account the nested structure of 
the A-EES, we propose reparametrizations of the type du* = g(u>i) dv that 
separate geometry from dynamic. Some of those parametrizations turn out 
to be generalization of the Jacobi amplitude. In Part II we consider geomet¬ 
ric properties of the A-system and the numeric computation of the functions 
involved. It will be published elsewhere. 

keywords: Integrable systems Generalized Euler system Jacobi and 
Weierstrass elliptic functions third Legendre elliptic integral 


1 Introduction 

We are interested in the real functions u>i(v) which are solutions of the integrable 
system of differential equations 

(1 <i,J<N), (1) 

with coefficients and initial conditions a,;,cnpO) el. Our study is based on the 
quadratic expressions 


Cij{y) = on LOj{y ) 2 - ajOJi{v ) 2 


( 2 ) 
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which are integrals of the system Initial conditions (IC) will be denoted 
cj° = w(0) = (wi(0),... ,cj n (0)). To simplify expressions we will use as notation 
Ui = uii(v) and = doji/dv. 

From the geometric point of view, the integrals ([ 2 ]) tell us that the flow defined 
b y 0 is the result of the intersection of quadrics in dimension N ; more precisely, 
elliptic and hyperbolic cylinders. Thus, the WEES family belongs to a larger 
family where the paraboloids are also included, as well as the degenerate cases 
defined by the hyperplanes. Its Poisson structure is defined by a determinant built 
on the gradients of the independent integrals, i.e. the Casimirs. When N = 3 
the classic mixed product is precisely the determinant: one of the integrals is the 
Casimir and the other the Hamiltonian; details will be given elsewhere [3]. 

One of the features of the system 0 is that it allows, from a dynamical sys¬ 
tem point of view, dealing with a large family of functions in the real domain 
in a unified way. It ranges from trigonometric functions (harmonic oscillator) to 
elliptic functions (pendulum and free rigid body), including also rational func¬ 
tions (for unbounded trajectories), etc. We will learn that different systems will 
allow us to introduce the same functions. For instance the hyperbolic functions 
may be introduced with N = 2, but also appear when N = 3 and two of the 
coefficients are equal). The interest of the study of the generic system N > 4 
(the case N = 4 is special, as we show below) lies in the fact that we face then 
hyperelliptic integrals and their inverses, a well established theory of special func¬ 
tions of complex variable made in XIX century which, nowadays, is in a revival in 
several branches of science, particularly in mechanics. But, although the theory 
is ‘at hand’, nevertheless its application results a nontrivial task, because of the 
number of parameters involved in the definition of the functions, solutions of an 
IVP. 

1.1 On Euler system, Jacobi functions and 3-EES 

In this paper, our program is to generalize Jacobi elliptic functions. Thus, within 
the dynamical system point of view we have adopted, let us remember how all 
this started. The history of the WEES begins with the well known Euler system 
of nonlinear differential equations in three dimensions [10] , giving the reduced 
dynamics of the free rigid body problem (the dynamics of the angular momentum 
vector II in the moving frame) 

ni = ai n 2 n 3 , n ' 2 = a 2 nin 3 , n ' 3 = 03 nin 2 , ( 3 ) 

such that ^2 oii = 0, where a* are functions of the principal moments of inertia. 

Associated with 0, the second fundamental system, known as the Jacobi 
system, is given by 

c n'i = 0 J 2 w 3 , J 2 = — W 1 W 3 , w 3 = -mwiw 2 , (4) 

with w(0) = (0,1,1). The functions solution of Q, denoted as aq = sn, w 2 = cn 
and cu 3 = dn, are called Jacobi elliptic functions. Then, the solution of ([ 3 ]) are 
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given by means of those functions, using the method of undetermined coefficients. 
For some readers couid be useful to consult our paper [2] where we have studied 
the extended Euler system 

= aiW2W3, oj' 2 = «2Wiw 3 , w 3 = a 3 a;ia;2, (5) 

i.e. the ([!]) for N = 3, considering generic values for coefficients ol{ and initial 
conditions defining the system. 

Relying on the work of Tricomi [13], Hille [7] and Meyer [IT] dedicated to sys¬ 
tem we have shown in a straightforward manner how Jacobi and Weierstrass 
elliptic functions in the real domain are connected with this system [2], although 
the tradition is to treat them separately because of the their intrinsic differences 
in the complex domain (see for instance Whittaker and Watson [14] and Lawden 
M)- Here we will apply the same approach to the system in ^dimensions. More 
precisely, we will present the generalization of both types of functions, where 
the N- Weierstrass function relates with the norm of the vector defined by the 
functions u;*. 

1.2 Integrals, functions and regularization 

Moreover, as an alternative to confront directly with hyperelliptic functions, we 
propose to experiment with reparametrizations starting from low dimensions. 
More precisely, we extend the regularization du* = w 3 du, already studied for the 
case IV = 3 by Molero et al. [12] , This way of proceeding seems to be an open 
line of work. The fact that elliptic and hyperelliptic functions are ‘naturally’ 
introduced within the context of complex functions may explain why we have not 
found references. It is due to the consideration of those functions in a dynamical 
systems context, in the real domain, that the regularization enters on the scene. 
More precisely we focus on ‘regularizations’ of the type dv* = g{uji)dv, a technique 
well known in classical fields such as Celestial Mechanics (where they are used 
for studies ranging from collisions to efficient numerical integration schemes). We 
will see that the new variable is a generalization of the Jacobi amplitude. This 
procedure, based on the symmetry of the system, alleviates the manipulation of 
the hyperelliptic functions involved, which are relegated to only one quadrature 
(the regularization equation ), separating it from the geometry (it is part of our 
research, knowing more on how generic this procedure is). 

This research has two parts. Part I, which makes the content of this paper, 
works in detail the cases N = 4, 5. The key aspect associated with this case is 
that for each IVP we deal with two or three parameters. In Section [2] we briefly 
refers to the equilibria as well as particular solutions such as the rectilinear. After 
that we fix the dimension considering the case 4-EES. In Section [3] we present a 
basic feature related to the ratios of the functions. In Section [4.21 we focus in a 
biparametric system, which we dubbed as Mahler system. In Section [5] we apply 
to our system the regularization technique. We identify that the new variable 
is a ‘generalized amplitude’. In Sect. [6] we provide with the addition formulas 
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associated to the Mahler system. Using them we propose extending the work of 
Bulirsch and Fukushima, we introduce some formulas related to the numerical 
evaluation of a 4-EES. In Section [7] we approach the system for N = 5, focusing in 
one of the particular cases, showing its connection with the previous dimension. 
Finally, as an application, we briefly consider in Section [9] the free rigid body 
formulated in Andoyer variables 

For the benefit of the reader we include two Appendices which contain prop¬ 
erties of 6i and elliptic Jacobi functions. There is a Part II, devoted to generic 
features of ([Tj) from the geometric point of view, and to the numeric evaluation 
of the Mahler system, following the steps of Bulirsch and Fukushima. This will 
be published elsewhere. 

We ought to close the Introduction pointing out that this paper does not 
contain a complete analysis of the relative role of the parameters involved in the 
defined functions. Some transformations related to the range of those parameters 
are required, similar to the well known transformations for the elliptic modulus 
of the Jacobi functions. That analysis is still in progress. 

2 Some basic features of iV-EES 

We have mentioned in the Introduction that our interest in this paper focuses on 
the study of some systems ([Tj) of low dimension. Nevertheless, as in any dimension 
common features are present, it is worth to briefly refer to some of them. 

2.1 On particular solutions: equilibria and straight lines through 
the origin 

Before we start our analysis of the IVP, a first question is to identify the equilibria 
of the system Q. Denoting P = (pi,P2, ■ ■ ■ ,Pn) an equilibrium point, we easily 
check that the system has the following set of equilibria: 

• Origin P = 0 E R n , 

• For n > 3, the points: Pi = (0,... ,pi ,..., 0), 1 < i < n, functions of the 

initial conditions. 

• For n > 4, planes n ilii2 = (0,... ,p ix ,... ,p i2 ,... 0), 1 < k < h < n, 

functions of the initial conditions. 

• For n > 5, the hyperplanes 

—2 (0) • ■ ■ ) Pi\ 1 • ■ • j Pi.2 1 • • • > Pin—2 1 • • • 0); 

1 <il<l 2 < in -2 < n. 

Thus, associated to these equilibria hyperplanes, we have the study of their in¬ 
variant manifolds and their connections, generalizing the heteroclinic trajectories 
in three dimensions. This is out of the scope of the present paper. 
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Straight-lines through the origin. Meanwhile in the generic study of the quadra¬ 
tures associated with our system (see Sect. 2.2) an assumption is commonly made, 
namely, the roots of the polynomials involved are different, when considering an 
IVP we may be under a scenario where we have multiple roots. This is precisely 
the case with straight-lines through the origin. Then, instead of requiring the use 
of special functions, the solutions are expressed by means of elementary functions, 
different for each dimension. 


2.2 Reduction to quadratures: Generalized Weierstrass function 

Taking into account the integrals and proceeding like in the classic case 
N = 3, we may reduce the system to a fundamental differential equation in two 
forms. The first one, after choosing one of te functions, say uji, it leads to the 
differential equation 


diUj \ 2 _ 3-N 

dv ’ i 


N 


[ )] . 


( 6 ) 


or, by separation, the corresponding quadrature 


(3—AO/2 
a) v = 


do 


n-M a ^ + c D} 1/2 


(7) 


As an alternative, if we introduce the square of the norm 

N 

n N (v) = uj(v) 2 = ^2uii(v) 2 , (8) 

i— 1 

after some straightforward computations we obtain 
i 0 2 n N 

(^= 4 n^-M, x>=°’ ( 9 ) 

2=1 2=1 

a differential equation whose solution Qn(v) may be seen as the generalized Weier¬ 
strass function p(v). Following either way we confront generically hyperelliptic 
integrals. 


2.3 On the normalized iV-EES 

Associated to a generic WEES 0, i.e. assuming that ^ a* / 0, we consider the 
square norm function ([8]) that satisfies 


dw 

dt> 


N N 

£>*)-IK 


2=1 2=1 


( 10 ) 
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Thus, introducing the functions 


Ui = 


OJi 
UJ ‘ 


we have 


which may be written also as 



d Cjj 
dv 


= c, 


N 

if 


.TV-4 


where the coefficients 


Ci = otiU 2 - (^ai)uj, 


( 11 ) 


( 12 ) 


(13) 


are integrals of the flow, whose values are determined for each IVP by the initial 
conditions. In other words, carrying out the reparametrization v —> v* given by 


d Gh 

dv* 


= c 


with initial conditions 


’ N ~ 4 dv, 

(14) 

i system 


N 



(15) 



W (°) 2 = 

(16) 


i.e. the flow (15) lives in and, like the differential system satisfied by the 

angular momentum in 3-D, we have ^ q = 0. Note that to deal with the system 
(15) versus (flj) will bring advantages, at least from the numerical point of view. 

With ( |15[ ) integrated we have = uj t (v*). Then, we still have to implement 
the quadrature associated to the regularization (14) in order to recover the rela¬ 
tion with the original variable. For instance, considering the first integral ci we 
obtain 


QU 


du* 


(17) 


whose quadrature gives the parametrization relation, solved generically by nu¬ 
meric methods. Note that, the case N = 4 is special, because we do not need to 
do regularization. 

Moreover, we will not pursue here with the study of the normalized system 
(15). For details on this analysis we refer to [3]. 

Let us close this Section pointing out another basic feature of this system; we 
refer to it as the scaling factor. If the functions u>i(v), (i = 1,... N) is a set of 
solutions, then taking a constant c, the functions Ui(v) = cu)i(c N ~ 2 v) satisfy the 
same system with the corresponding IC given by m( 0) = cwj(0). We will make 
use of this property along the paper. 
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3 The case N = 4. Relying on Jacobi elliptic func¬ 
tions? 

We focus now on the 4-EES case. For each IVP, with some abuse of notation, 
we refer to the functions solutions generically with Wj. Later, referring to some 
specific systems, we will introduce new notations. 

At this point, perhaps some readers would like to know the original motivation 
of our interest in 4-EES case. The reason is connected with an observation about 
the classical way in which the study of the rigid body dynamics is developed, 
based on Jacobi elliptic functions. Meanwhile those functions depend on one 
parameter (elliptic modulus), and appear naturally tied to problems like the 
pendulum or the measure of an arc of ellipse, when we apply them to the rigid 
body problem, we need to consider a second parameter (the characteristic, a 
function of the principal moments of inertia). In other words, the first and third 
Legendre elliptic integrals are involved. Since Jacobi, the way to proceed has 
been: (i) to introduce complementary functions Z and 0; (ii) to make use of the 
addition formulas of elliptic functions, dealing with the second parameter as an 
amplitude, etc. Here we search for an alternative to such approach considering a 
generalization of Jacobi elliptic functions with two parameters. 

Thus, we start with the 4-EES 

w'l = aiCU 2 W3W4, 
u' 2 = a 2 ijJi W3W4, 

CJ3 = 0:3 CJl W2W4, 

W4 = a 4 o;i w 2 w 3 , 

with given initial conditions <u°, and the corresponding six quadratic first integrals 
Q, of which three are independent (Fig. [I] shows a graph of the solution of the 
system ©)• Although by scaling and a change of variables we could get rid of 
two of the coefficients cq, for our purpose it is convenient here to maintain all of 
them. 


(18) 



To our surprise, the only reference we have found so far to ( |18| ) is E. Hille [7], 
where the case N = 4 is considered in Chapter 2 (exercises 7, 8 and 9) under the 
suggestion of K. Mahler. More precisely he considers the IVP cu(0) = (0,1,1,1) 
and coefficients a.% = (1, —1, —oi 2 , — /3 2 ), with both coefficients less than one. He 
says “the solutions are hyperelliptic functions of genus 2”, a statement on which 
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we disagree. Finally he mentions “the example can be generalized in an obvious 
manner.” 

Thus, our plan is: (i) to study (18) as an extension of the case N = 3 where the 
Jacobi elliptic functions were defined. Note that represent a drastic reduction in 
the number of parameters (coefficients and IC) to discuss; (ii) To introduce again 
regularizations. In order to approach both aspects, apart from its own interest, 
we think the case IV = 4 is critical in the search for methodologies to follow 
when dealing with systems of higher dimension, i.e. in the reign of hyperelliptic 
integrals. 


3.1 Nested structure and integration by Jacobi elliptic functions 

Extending what we know for the case N = 3, a basic feature of the 1V-EES is 
its relation with the system verified by the ratios. Referring to that we say the 
TEES has a ‘nested structure’, and we call it the ‘Glashier Ratios Property’. 
Moreover the case N = 4 asks for a particular study devoted to it. As we will 
see, for other dimensions a regularization is needed. 


Proposition 3.1 (Glashier Ratios Property) Given the functions ufiv) ver¬ 
ifying a 4-EES, then the functions Ui(v)/uj(v) defined by their ratios, ( i,j , k, l ) E 
Per(l,2, 3,4) satisfy a 3-EES given by 


d_/(T\ _ c i Wfc 
d v \uiJ 1 ui wfi 
d_/wA = c i 
d v \LOiJ 3 uy ufi 

d /cu fc \ = C ! UJj, UJj_ 

d v \ui J k uii u>i ’ 


(19) 


with initial conditions wfiO) / wfiO) and coefficients given by the integrals C\ = 
ctiuf - aiwf. 

Proof : It is straightforward making use of the definition of the 4-EES q.e.d. 


Remark 3.1 From the previous Proposition \3.1\ readers familiar with the expres¬ 
sions of Jacobi elliptic functions, and their computation by means of Jacobi theta 
functions 9 fix), may wonder what the relation between those functions and the 
ufiv) might be. We have gathered some of those systems in an Appendix. In fact 
the reader will find in Lawden (Chp 1) a number of properties ofQi functions which 
are also satisfied by the w*. Perhaps, the simple fact that 0\ (0) = 02(0)03(0)04(0) 
is satisfied for the 4~PPS when we take a\ = 1, is one of the most surprising. 
We will come back to this below. 

Remark 3.2 Note that there is the possibility to take a slight different version of 
the ratios, namely to work with u* = CjUi/ojj, with coefficients c) still to be de¬ 
termined, in order to simplify some expressions, adjust constants in applications, 
etc. We do not follow this alternative in this paper. 
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Proposition 3.2 For suitable IC the 4--EES (18) has as solution the bounded 
functions uiiiy) = 0 )) given by 


wiO) = Cf 
u 2 {v) = Cf 
W3(u) = Cl 


sn(av\m,i) 

\Jl — n\ sn 2 (av|mi) ’ 
cn(av|mi) 

V* — n± sn 2 (av|mi) ’ 
dn(au|mi) 

A — 77,1 sn 2 (a?;|r 77 ,i) ’ 


m(v) = Cf 


A — ?7 1 sn 2 (a?;|777,i) ’ 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 


where sn(au| 777 i), etc are the Jacobi elliptic functions, and the constants Cf, a, 
777i and 771 are functions of a* and 0). 

Proof : Let us assume IC cu° = (a;) 1 ,... , 074 ) such that Uj 7 ^ 0 in its domain of 
definition. According to the previous Proposition, we consider the ratios and the 
reciprocals 1 /uj, that we denote 


j ■ _l ■ 

K = ~, 


UJj 


j 

u j = 


(24) 


in the domain where ujj is defined. Without loss of generality we assume we refer 
to the case j = 4, with IC such that 104 > 0. Moreover, we still simplify a bit 
more the notation writing uf = U{. 

Then, according to Proposition 3T it results for the functions Ui , i = 1,2,3 
we have the following system 


u\ = Cfu 2 u 3 , 

u ' 2 = ClusUi, (25) 

U 3 = C 3 U\ U 2 , 

with IC Ui( 0) = u® = 0 J®/ujj . Moreover, from the first integral 

ot\w\ — 0 . 41^1 = Cf (26) 


we may write 


u l = -a 4“i)- 


(27) 


Because the functions Ui, i = 1,2,3 satisfy (25), they belong to the set of func¬ 
tions defined by the ‘Jacobi elliptic functions’ sn, cn, dn and their ratios. Then, 
following Crespo and Ferrer [2], we know our system corresponds to one of the 
four possible cases (Glashier systems), depending on the sign of the integrals. 
Here, to continue our reasoning on the system (18), we focus on the case where 
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the sign of Cf is different of C 2 and C| (the other cases are treated likewise). 
This means that Ui , i = 1, 2, 3 are of the form, say 


u±(v) = <5i sn(au, mi), 

u 2 (v) = S 2 cn(av,mi), (28) 

u 3 (v) = f ?3 dn(au, rn\). 


Proceeding by the method of undetermined coefficients, replacing (28) in (25) 
we identify that the constants Si, a y mi satisfy a system of algebraic equations 
whose solution is 


S 2 = u° 2 , S 3 = v%, <5i = \J-oli /a 2 S 2 , 

a = ai5 2 5 3 /5i, mi = ct 3 Sj/(a 2 Sl) 

(for details see for instance Lawden (9] , p. 132). 


Summarizing, according to 


and 


we have uoi = m/u 4, where Ui 


(i= 1,2,3) are the Jacobi elliptic functions and U 4 is given by (27). From those 
expressions, we obtain the functions p0l)-(23), where 


(29) 


C\ = ^Cf/a 1 , Cf = Si/Cj, m = a^l/ai 


and, as stated in the Proposition, initial conditions still have to be chosen such 
that ni < 1- q.e.d. Before we continue it is convenient to 

formulate the previous Proposition in a ‘complementary form’, where we make 
more transparent the role played by coefficients and initial conditions. 

Proposition 3.3 The functions uoi(v), i = 1, ... 4 , given by 

w 2 ( 0 ) w 3 ( 0 ) u> 4 ( 0 ) sn(au|mi) 


oji(v) = 
oj 2 {v) = w 2 (0) 
iv 3 (v) = w 3 ( 0 ) 
U)i(v) = CJ4(0) 


o y/l + ni sn 2 (at>|mi) 

cn(av\mi) 


y/l + ni sn 2 (av|mi) 
dn(ov|mi) 

i/l + ni sn 2 (at>|mi) 
+ ni sn 2 (au|mi) 


satisfy a differential system of the type (18) given by 


(30) 


UJ 1 = UJ 2 C 0 3 CO4 
U }' 2 = -(1 + ni ) 


^(0M(0) 

/ / s a 2 

co 3 = —(mi + n 1 ) 


'^(0)a;|(0) 


wi W3W4, 


CO 1 C 0 2 LO4, 


004 = -n 1 - 


■Wl co 2 oo 3 , 


with co = ( 0 ,cj 2 (0), w 3 ( 0 ), ^ 4 ( 0 )) as initial conditions 


( 31 ) 


10 






















Proof : It is a straightforward exercise by computing derivatives. 


q.e.d. 


Remark 3.3 In particular, choosing uii( 0) = 1 (i = 2,3,4) and a = 1, join with 
n\ = n and = m — n in Proposition\3.3[ we have the Jacobi elliptic functions 


sn(u) = 


ui(v) 


cn(u) = 


u- 2 (v) 


dn(u) = 


w 3 (u) 


CJ4(u)’ V ’ UJ^v)' V ' UJ/±{v) 

with elliptic modulus m\ = m — n, where w,(r: m, n ) satisfy the system 

Oj'l = W2W3CJ4, 

uj' 2 = —(1 + n) UJ3 W4, 

u>' 3 = —m coi W 4 , 

W4 = —nuj\ W2W3, 


(32) 


with integrals 


uj ‘2 + (1 + n) uj\ — 1 , 

o o 

UJq + TYIUJ-^ — 1 , 

W4 + nuj\ = 1. 


(33) 


If 0 < n < m < 1, we have —1/^/1 + n < uj\ < l/\/l + n, — 1 < W 2 < 1, 

— m /(1 + n) < CJ 3 < 1 and \J\ — n/(l + n) < W 4 < 1 . 

More details on the system (32) will not be given in the rest of this paper. 


4 Studying two 4-EES systems 

Looking for the generalization of Jacobi elliptic functions, we now focus on two 


cases of (18): 


One-parameter (#j similar) family in Sec. |4.1| and 


Two-parameter family (Mahler system) in Sec. 4.2 


It is worth noting that the first two equations in both systems (see (38) and (39)) 
are equal, with the consequence that one of the integrals is uf+uj^ = 1, which is 


not the case for the previous system (32). 


In relation with both, before we continue, a comment on notation is due. In 
what follows, it is convenient to redefine some of the constants which appear in 


the previous expressions. More precisely, in Sec. 4.1 we write mi = k 2 , and we 


will find that a and ni are functions of k. Likewise, in Sec. 4.2 we fix all initial 


conditions and coefficients except two of them, denoted by —m and —n. 
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4.1 One-parameter uJi(v) functions, ‘similar’ to Jacobi 6i functions 


We look here for functions Wj, solutions of our differential system (18), similar 
to Jacobi 9i functions. What we mean by that should be made more precise: (i) 
coefficients and IC of the 4-EES have to be dependent only of one parameter: 
oti = ati(k), = w°(/s); (ii) Moreover those functions u n(w,k) ought to be found 

imposing that they verify properties defining 9, de Jacobi. 

Such search does not appear straightforward because, we remember, 9i func¬ 
tions are defined as 1-parameter Fourier series solving the heat equation. Our 
way of proceeding will be to take into account those properties of 9i which could 
be imposed on the differential system: both the ratios and the identities satisfied 
by #j(0) are essential for us. 


Proposition 4.1 (w*: ‘similar Jacobi 9% functions’) Choosing initial condi¬ 
tions as functions of the elliptic modulus 


wi(0) = 0, uj’2 (0) = Vak , uj 3 (0) = \Ja, u; 4 (0) = V a k' 


join with 


2 K 


a =-, n\ = k' — 1, 


7 r 


m\ = k 2 


where k' = \/l — k 2 , then we may write 


ui(wjj(0)z) 

v 2 (uj%( 0 )z) 

W3(wf(0)z) 


^3(0) ui(z) 
o; 2 (0) w 4 (z)’ 
^ 4 ( 0 ) U2{z) 
o; 2 (0) u 4 (z)’ 
^ 4 ( 0 ) U 3 {z) 
u; 3 ( 0 ) oj 4 (z) 


(34) 

(35) 


(36) 


in other words, we express the Jacobi elliptic functions as ratios of the w«(u), in 
a similar way as Jacobi gave them with respect to the 9i functions. 

Proof : .- It is a straightforward exercise replacing the previous values (34) and 
(35) in Proposition 3.3 The result is that the functions are 


lo\(z, k) 
td 2 (z, k) 
w 3 (z,k) 
u 4 (z, k ) 


snm 


Va kk' —7=_, 

\J\ — (1 — k') sn 2 (u) 

c "(“> 

s/\ — (l — k') sn 2 (u) 
dn(rt) 

\Jl — (1 — k') sn 2 {u) ’ 
yjl — (1 — k') sn 2 {u) ’ 


( 37 ) 


join with u = az. 
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Thus the system (31) given by 


uj 1 = UJ 2 CO 3 U 4 , 


u 2 = 


u 3 = 


<^4 = 


-w 3 w 4 Wl, 

1 - k' 


k 

l-k' 
~k 


■ cu 4 u 1 w 2 , 


-Wl CJ 2 w 3 , 


( 38 ) 


with initial conditions (34), is the IVP we were looking for. Fig. [ 2 ] shows an 
example of a graph of this set of functions. q.e.d. 
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Figure 2: Graph of the ^-similar for m = 0.95. 


It is an exercise to check that the functions (37) verify identical relations to 
the linear combinations satisfied by the square of Jacobi (9,; functions (see Lawden, 
formulae (1.4.49)-(1.4.52), p. 11). 


4.2 Mahler system. A biparametric 4-EES: 

As a second distinguished 4-EES we consider now a ‘biparametric’ case we call 
Mahler system. It is an IVP which defines the functions uji(v\m,n), solutions of 
(18) depending on two parameters, such that 

• coefficients a = (1, — 1, — m, — n ) 

• initial conditions uj° = (0,1,1,1). 


When n = 0 then Ui(y) are the Jacobi elliptic functions and = 1. 

Note that this represents some abuse of notation, because n has already been 
used to denote the last component of an A-dimension system. Nevertheless, we 
think by the context it will become clear when is a coefficient: n G M, although 
in some occasions n might be used as a counter (ordinal number: n £ N). 


Proposition 4.2 (Mahler system) 

The 4-EES given by 

0j[ = UJ2UJ3UJ4, 

U}' 2 = —Ul CO3 OJ4, 
uj' 3 = -mwi w 2 0J4, 
W4 = —n U>1 LO2LO3, 


(39) 
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where n < m < 1, with IC w(0) = (0,1,1,1), has the functions 

sn(av\mi) 


u>i = A 


U! 2 = 


\J\ — m sn 2 (au|mi) ’ 
cn(au|mi) 


w 3 = 


UJ 4 = 


0 “ — n\ sn 2 (au|mi) 

dn(au|mi) 

V 7 ! — n\ sn 2 (au|mi) 


( 40 ) 


0 " — n\ sn 2 (av\mi ) ’ 
as solution, with values a, A,m±,n± given by 

a = \/l — n, A = 1/Vl — n, 


n 


n l = 


n — 1 


mi = 


n — m 
n — 1 


(41) 


Proof: Let us consider the system (39) as an IVP with lo (0) = (0, u; 2 (0), w 3 (0), W4(0)), 
(wj(0) / 0, i = 2, 3, 4) dependent of two parameters (m, n). It admits as solution 
the functions 


u)i(u) = A 


sn(au|mi) 


w 2 (u) = 07 2 (0) 
Cb 3 (v) = w 3 (0) 
£ja(v) = OJ4 ( 0 ) 


\fl — ni sn 2 (au|mi) ’ 
cn(an|mi) 


0 " — ni sn 2 (an|m-i) ’ 
dn(an|mi) 

v/i" — ni sn 2 (an|m.i)' 
1 

v/i" — ni sn 2 (an|mi) ’ 


where a, A, mi, ni are given by 


a = uj 3 (0 )^ujI(0 ) -nw|(0), 
w 2 (0) o; 4 (0) 


A = 

ni = 

mi = 


y/uj(0) - nw|(0) ’ 
nu |(0) 

nuj |(0) — w|(0) ’ 

w l(0) ( nw 3 (0) — mcnKO)) 


(42) 


w 3 (0) (ncu|(0)-u;f(0)) 
and the derivatives at the origin satisfy 

wi(0)' = w 2 (0) w 3 (0) cn 4 (0), Wj(0)' = 0, 


(43) 


where i = 2,3,4. Then, choosing as IC the quantities w(0) = (0,1,1,1) and 
replacing them in (42), we readily obtain the values © for those parameters. 

q.e.d. 
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Remark 4.1 In particular the case n = 0 leads to: a = 1, A = 1, n\ = 0 y 
mi = m, i.e., the Jacobi elliptic functions. We have another special case when 
m = 0. As we have assumed n < m, in this case n < 0 and the diferencial system 


(39) corresponds again to a Jacobi system, but now with negative parameter (there 


is a transformation to reduce it to the normal case, see Appendix B, Sect. 0 - 
For more on particular cases see Section 5.f. We leave for the reader to work 


out the other particular cases defined by special values of the pair (m, n). 


5 Regularization and ‘generalized amplitudes’ for the 
Mahler system 

We have just solved the system N = 4 in the standard way: making use of known 
functions (Jacobi elliptic functions). In what follows we are going to proceed 
making use of the regularization. To do that, we start remembering in Section 
5.1 the recent proposal of the authors for N = 3 (see Molero et ol. [12] 1. which 


is intrinsically connected with the Jacobi amplitude. After that we develop the 
same approach for the N = 4 case. That proposal entails to study, at least, two 
possible regularizations v —> v* given by 

• dv*/dv = W 4 , 

• dv*/dv = CU 3 OJ 4 , 


which we gather in Sections 5J2 and 5JI Let us proceed one by one. But, before, 
we remember in Sect. 15.11 how this has been done for the 3-EES. 


5.1 Preliminaries: 3-EES and regularization 

Let us consider the 3-EES |5]) with initial conditions cv° = u( 0) = (cji(0), W 2 ( 0 ), <^3 (0)), 
whose values we choose below. This system has the integrals 


O 0 0 O OO 

a\tjj 2 — *^2 = Cj, 01^3 — 0:3 u>i = C\ . 


(44) 


Let us assume ccj and IC such that uis(v) > 0. Then, making use of the 
parametrization 

rh,* 

(45) 


the system ([ 5 ]) reduces to 


dwi 

dib 


du 


— «1 


= w 3 , 


di02 


du* 


— a 2 wi, 


(46) 


join with the quadrature defined by (45). Choosing the coefficients a\ = 1 


«2 = -1 and IC (wi(0), 0 ^ 2 ( 0 )) = (0,1), the system (46) defines the trigonometric 
(circular) functions: 

sin(u*), cos(w*). (47) 
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(with other conditions, by a change of variables we may reduce it to this case) 
Then, keeping in mind (44), the regularization (45) takes the form 


dv* 

dw 


— \J ' C i + CK3 w l- 


(48) 


Motivated by the dynamical system defining the simple pendulun ^ it is chosen 
CJ 3 (0) = 1 join with 03 = — k 1 2 , where k 2 < 1. Thus, replacing in (48) we have 


dv = 


dv* 


\Jl — k 2 sin 2 v* 

whose quadrature and inversion leads us to the Jacobi “am” function: 


(49) 


v* = am(u, k). 


Finally, replacing in (47) we have the Jacobi functions 


sin(u*(v)) = sin(am(i>, k)), 
cos(v*(v)) = cos(am (v,k)), 


which today, following Gudermann, are denoted in the form 

sn(v; k) = sin(am(v, k)), 
cn(u; k) = cos(am(u, k)). 


(50) 


(51) 


Completing our set of functions W 3 is given by 

W 3 (v) = dn(u; k) = \/l — k 2 sn 2 (u; k). (52) 

Summarizing, using the previous notation, the integrals ( |44[ ) lead us to the well 
known expressions relating these functions 

sn 2 + cn 2 = 1, dn 2 + fc 2 sn 2 = 1. (53) 

Finally, replacing in ([ 5 ]) we write what some authors refer as “derivation rules ” 
of Jacobi functions: 


sn 7 = cndn, cn 7 = —sndn, dn 7 = —fc 2 sncn. 


(54) 


5.2 The dv*/dv = (u 4 regularization. 

Proceeding as in the previous Section we treat now the case IV = 4 by means of 
the regularization 

dV * /rro 

= uj 4 . (55) 

du 


1 This lead us to an interpretation of the regularization: v = t and v* = 4>, in other words 

‘time’ and ‘angle’. Angle in the 1-2 plane; arc through the integral uif + = 1, a circle 

projection of the integral which is a cylinder. 
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Remark 5.1 Remember the comment above in relation with notation; although 
there is some abuse using again v* for denoting the new independent parameter, 
from the context we distinguish it from the one studied in the previous Section. 


As a consequence the system (18) is reduced to 


duq duq duq 

— = 01^2^3, — = «2WlW3, —— =asUi(W2, 

dv* du* dv* 


and uq(u*) which will be obtained using one of the integrals, after we have solved 
the previous system. 

We focus on the case aq = 1, aq = — 1 and = —m because, as we have 
said before, we plan to generalize Jacobi elliptic functions. Thus, we have 


uq = sn(t;*;mi), uq = cn(u*; mi), uq = dn(w*; mi) (56) 


and for the differential relation using the integral n uf + = Cf and the initial 

conditions, we may write 


v = 


dv* 

y/l — n\ sn 2 (u*; m-i) 


(57) 


5.3 IV = 4. The regularization dv*/dv — uq 0 J 4 . 


Proceeding the same way as for N = 3, we treat now the case N = 4 by means 
of the regularization 

du* 

=w 3 w 4 . (58) 

du 


As a consequence the system (18) reduces to 

doq 


du* 


— Ot 1 W2, 


dcj 2 

du* 


— a 2 wi, 


(59) 


and two quadratures associated to W 3 y W 4 . In fact, they are not needed because 
the integrals gives us 

Jj = C[ - oqwi, (i = 3,4). 

Note that C\ are constants which depend on the initial conditions. 

Without loss of generality we will assume our system is made of bounded 
functions. Then, by a change of variables, our system (59) reduces to aq = 
1 , «2 = — 1 , thus it results 


cui(u*) = sinu*, uj 2 (v*) = cosv*, 


Considering the previous integrals we may write (58) as follows 

, du* 

dn = 


^Uc\ 


(60) 


( 61 ) 


— cq sin v 
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or in a slightly different form 


Adw 


du* 

y/{\ — /3\ sin 2 u*)(l — sin 2 v*) 


(62) 


where f3i and A are functions of C\ and a*. 

In what follows, with the Mahler system in mind as the basic 4-EES, it is 
convenient to take the associated notation: 


j3i = n, 02 = m, A = 1. 


In other words, the differential relation (62) reads 
du = 


du* 


\/(l — n sin 2 u*)(l — msin 2 v*) 
The quadrature takes the form 

dd 


v = 


r v 

G(v*,n, m) = 

Jo 


!o \/(1 — n sin 2 i?)(l — m sin 2 d) 
Thus, we dehne the period as the two-parameters function 

T / 2 dd 


G(ir/2, n, m) = 


'o \/(l — nsin 2 i?)(l — msin 2 19) 


(63) 


(64) 


(65) 


Thus, when (n, m) = (0,0), we have G( 0,0) = 7t/ 2, and when (n,m) = (1,1), we 
have G(l, 1) = oo. 

When (m, n) are small, if we carry out the Taylor expansion of the integrand, 
after the evaluation of the quadratures, G(n, m) may be approximated in the 
form 


G(n , m) = 

7 r r, m 9 m 2 25 m 3 1225m 4 

— 1 T — T-T-T- 

2 L 4 64 256 16384 

n/ 3 m 15 m 2 175 m 3 2205m 4 \ 

+ 4 V 1 + “8~ + 64 + 1024 + 16384 ) 

9n 2 / 5m 35m 2 105m 3 2695m 4 \ 

l 1 + T2~ + 128 + 512 + 16384 ) 

25n 3 / 7m 189m 2 231m 3 3003m 4 \ 

+ 256 V + 16" + 640 + 1024 + 16384 / 

1225n 4 / 9m 99m 2 429m 3 6435m 4 y 

+ 16384 V 1 + ^0 + 320 + 1792 + 32768 /. 
+h.o.t. 


although the previous expression may be written in different form making more 
explicit its symmetric character with respect to m and n. 
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Now we define the generalized amplitud amg as the inverse function 

v* = amg(w; n, m). 


Thus, considering the expressions (60), we have 

sin v* = sin amg(r, n, m) = sng(v, n, m) 


( 66 ) 

(67) 


and 

cos v* = cos amg(r, n, m) = cng(u, n, m) (68) 

• There is an alternative way of proceeding. If we consider the change of variable 
sint? = x it allows to follow the steps of Jacobi for the case N = 3. Then, the 
differential relation (63) takes the form 


du = 


dx 


y/(l — x 2 ){\ — nx 2 )( 1 — mx 2 ) 


or, inverting the expression 

^ = y/(l — x 2 )(l — nx 2 )( 1 — mx 2 ). 

In other words, we define the function sng 

x = x(v, n, m) = sng(u; n, m ) 


(69) 

(70) 

(71) 


as the two-parameters function (whose range is made more precise below), solu¬ 
tion of the differential equation 

J = (1 — x 2 )(l — nx 2 )(l — mx 2 ). (72) 

In this paper we will restrict to a range n < m < 1. 

• Then, associated with sng we propose the following functions 

cng(u; n, m) = ±\A ~ sng 2 (u; n, m), (73) 

dng(u; n,m) = y /1 — m sng 2 (u; n, m), (74) 

fng (v;n,m) = \Jl — nsng 2 (t>; n, m). (75) 


To simplify notation we will write sng(u; n, m) = sng, etc. Examples of the graph 
of these new functions can be seen in Figs. [3]and[4j 

Due to the process we have followed, we immediately check that these func¬ 
tions sng, etc verify the following I VP 


dsng 

du 

dcng 

du 

d dng 

dv 

dfng 

du 


eng dng fng, 
—sng dng fng, 

—m sng eng fng, 
—n sng eng dng, 


(76) 
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with initial conditions (0,1,1,1). The integrals, as we have mentioned before, 
lead to the following expressions 

00 -iO O o O /\ 

eng + sng = 1, dng + msng = 1, fng + nsng = 1. (77) 


Thus, from the functions solution of the Mahler system, the Jacobi functions are 
given by 


sn(au; m\) 
cn(au; mi) 
dn(au; mi) 


1 sng(u; m, n) 
A fng(u; m, n ) ’ 
cng(u; 7n, n) 
fng(u; m , ?i) ’ 
dng(u; m, n) 
fng(u; m, n) ’ 


(78) 


• Taylor expansions of sng, eng, dng and fng near the origin. 


As a direct application of the definition of those functions by the differential 
system (76), we may easily compute to any order the Taylor expansion of the 
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previous functions: 


. . 1 + m + n o 

sng(u) = v - v 

6 

1 + 14 (m + n + mn) + m 2 + n 2 ^ 

H----- v b + ... 

120 

, s 1 2 1 + 4m + 4n 4 

cng(u) = 1 - - v +- — - v 

1 + 44(m + n) + 16m 2 + 104mn + 16n 2 6 

720 V + "' 

, . . m 9 m(4 + m + 4n) d 

dng(u) = 1 - - v 2 + —^-—-- u 4 

m(16 + 44m + m 2 + 104n + 44mn + 16n 2 ) 

- v + 

720 

n 2 , n(4 + u + 4m) 4 
fng(u) = 1 - - u 2 + ^ Z u 4 

n(16 + 44re + n 2 + 104 m + 44mn + 16m 2 ) 6 

720 W + ' 


( 79 ) 


Remark 5.2 The interest of these expansions is connected with the computa¬ 
tion of these functions. By extension of the process followed by Bulirsch and 
Fukushima computing Jacobi elliptic functions (see Appendix). Nevertheless, 
there is still work to be done comparing that scheme with the possible advantages 
of using regularization. 


5.4 Particular cases 

• ra = 0. In this case, due to the choice of the initial conditions, we have 
fng(u) = 1. Moreover we have sng(u;0, m) = sn(u,m), etc, i.e. the Jacobi 
elliptic functions with elliptic modulus m. 


m = 0. Here, based on the initial conditions, we have dng(u) = 1. Moreover 
sng(u; n, 0) = sn(u, n), etc, i.e. the Jacobi elliptic functions have an elliptic 
modulus n (que es negativo, thus we still needs to make a transformation; 
see (41) leading to mi). 

m = 1. In this case the differential equation is 


— - = (1 — x 2 )\/l — nx 2 . 
dv 

For this quadrature we obtain 

_ 1 (1 + x) (1 - nx + a/(1 - n)( 1 — nx 2 )) 

2^/1 — n (1 — x) (1 + nx + a/(1 — n)(l — nx 2 )) ’ 

whose inversion is possible, because it is injective. 


(80) 


( 81 ) 
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• m = n. Now the differential equation is 


= (1 — mx 2 )\/1 — x 2 . 
av 


We obtain 


v = 




m 


-ArcTanl vl — m 


Vl - • 

Again the inversion is possible because it is injective 


tan(\/l — mv) = \/l — 


m 


VT^. 


More precisely, we have 


(82) 

(83) 


(84) 


x = 


tan(v / l — mw) 


sjl-m + tan 2 (Vl — mv) 
Graphical examples for n = m can be seen in Figs. [5] and [6j 


(85) 



1.0 

0.8 

0.6 

0.4 

0.2 



V 


Figure 6: Mahler m = n = 0.95. Falta otra con valor mas extremo de n = m 


• m = n = 1. In this case 

and finally, after inversion, it results 

v 

x = , 

Vl + v 2 

• m = n = 0. In this case we recover the circular functions. 

There are other particular cases related to unbounded trajectories, like the 
straight-lines which are expressed by elementary functions. This requires 
the signs of the coefficients to be the same, something that we have excluded 
when choosing our system. 


( 86 ) 

(87) 
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6 Addition formulas 

In order to alleviate the notation, we introduce the following convention 

sng(ax; m, n) = s^, cng(ax; m, n) = c ax , 
dng(ax; nr, n) = d ax , fng(ax; m, n) = f^. 


Theorem 6.1 (Addition-Subtraction formulae for the 4-Mahler functions) 

The addition and subtraction formulae for the 4-Mahler functions are given next. 


sng(x ± y; m, n) = 
A (s, 


( 88 ) 


d aV f aX ± S ; 


ax ^ay ^ay X ax 


ay'-'ax 


dax fay) 


\J (faxfay m l^ax^ay)^ n l (s a xC a yd a yfax i S a yC a xd a xff 


ayy 


cng(x d= y; m, n) = 



Cay fax fay ~F S; 


^ax'-'ay 


^ax^aydaxday 


^•l^ax^ay)^ ^1 (SaxC a yd a yfax i S a yC a xd a xfay)^ 


dng(x d= y; m, n) = 

daxdayfaxfay ~F Sax^ayCax^ay 

\J (faxfay m l^ax®ay n l (SaxC a yd a yfax i SayCaxdaxfay) 


fng(x ± y; m, n) = 


f2 f2 ,-q g2 „2 

ax ay 1 °a.x °ay 



where A, a, m\ and n\ are given in formula (43) (en la proposicin 5). 

Proof : Let us prove the formula corresponding to sng(x±y; m, n), the remaining 
ones are analogous. By Proposition 5 we have that 

sn(ax + ay ; mi 
— nisn 2 (ax + ay ; mi) 

Thus, using the addition and subtraction formulae for the Jacobi elliptic sine (see 
Appendix B) and assuming the following convention 


sng(x ± y; m, n) = A 


sn(ax;mi) = s x , cn(ax;mi) = c x , dn(ax; mi) = d x , 


we obtain 


sng(x ± y; m, n) = 


A- 


SxCydy lit SyC x d x 
1 — mis 2 s 2 


'(1 - mis^Sy) 2 - ni(s x c y dy ± s y c x d x ) z 


(1 - ?nis 2 s 2 ) 2 
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simplifying denominators 

sng(x ± y; m, n) = 

A 


SxCydy i SyCxdx 


(1 - mis|s 2 ) 2 - ni (s x Cydy ± s y c x d x ) 2 


(89) 


Finally, recalling that 


s x = 


Cx — 


d x = 


1 sng(ax; m, n) 
A fng(ax; m, n) 
cng(ax; m, n) 
fng(ax; m, n) 
sng(ax; m, n) 
fng(ax; m, n) ’ 


and likewise for s y ,c y ,dy, if we multiply numerator and denominator in (89) by 
fng 2 (ax;m, n) and fng 2 (ay;m, n) we obtain (88) after algebraic simplifications. 

q.e.d. 

Corollary 6.2 The formulae for the double angle of the 4-Mahler functions are 
given by 

2 A S a xC a xdaxfax 


sng(2x; m, n) = 
cng(2x; m, n) = 
dng(2x; m, n) = 
fng(2x; m, n) = 


]/(& ~ m l s tx ) 2 - n l ( 2 SaxCaxdaxfax)' 

„2 r2 -j- „2 j2 
^ax ax i °ax u ax 


]/ (f^c - nil sf x ) 2 - 111 (2 SaxCaxdaxfax)' 

j2 c2 -j- g2 „2 

u ax i ax i °ax u ax 


(90) 


\J (f^x - m l s ax) 2 - n l (2 SaxCaxdaxfax)' 


r4 4 

f ax - m l s ax 


yf (f^x - m l S U 2 - 11 1 ( 2 SaxCaxdaxfax)' 


Corollary 6.3 The formulae for the half angle of the 4-Mahler functions are 
given by 


,x 


sng( —; m, n) = A 


fax C ax 


fax d ax Hi (fax Cax) 


,X 


cng(-;m,n) = 


dax + C ax 


,X 


dng(-;m,n) = 


fax d ax Hi (fax Cax) 

' (Cax + d ax )(fax + d ax) 

(fax + C ax ) (fax + dax) - ni(f 2 x - C 2 X ) 


(91) 


,x 


fng(-;m,n) = 


fax d ax 


fax d ax H].(f ax C ax ) 
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6.1 On the numerical computation of cUj functions by extending 
Bulirsch-Fukushima method 


As we know Jacobi elliptic functions are defined by some ratios of 9 t Jacobi 
functions. This way of handling the Jacobi elliptic functions is convenient due 
to the fast convergency of those series. Nevertheless, at present, fast numeric 
codes compete with this classical analytic approach. More precisely, in order to 
implement those codes addition formulas compute Jacobi elliptic functions are 
basic expressions in that process (see FukushimaO E]). 

We can extend those expressions to the cj, functions. Thus, as Fukushima 
explains, the algorithm is made of three steps: 

(i) the forward transformation defined by (Corollary 2, Half arguments formulas: 
© reducing the values of cj* by a number of iterations; 

(ii) evaluation of the Mac-Laurin series expansions given by (79) and; 


(iii) the backward transformation (Corollary 1: Double arguments formulas (90)) 
as many times as the forward transformation. 

Details of the implementation of this process will be given in |3j. 


7 On the case N = 5 

As we have pointed out in the Introduction, hyperelliptic integrals appear in Q 
when N > 5. Thus it is convenient to see in some detail the case N = 5, the 
lower system belonging to this category. 

Thus, as before, we start keeping the notation used in lower dimension 

= ai(j02CU3U4UJ5, 

J 2 = a 2 u 1 W 3 CU 4 W 5 , 

CU 3 = «3 0 J\ W 2 W 4 CJ 5 , (92) 

CO 4 = a 4 u;i w 2 co 3 u> 5 , 
u ' 5 = a 5 u>i W2W3W4, 

with given initial conditions w(0). As examples in Figs. [7] and [ 8 ] we present two 
set of functions of the 5-EES family. 


Uj 


1.0 





- COS - 


0.5 


— <U| 

^3 




-0.5 


2 


\ 6 / 

8 

/TO 









Figure 7: 5-Mahler system graphs for p = 0 . 2 , n = 0.4, m = 0.7. 

We will proceed as in the lower dimensions N = 3,4, considering alternative 
procedures to the classic solution based on the direct reduction to hyperelliptic 
integrals. In other words: 
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(i) we introduce the functions vPAv), (where we maintain the notation) ratios of 
the oji 

( 93 ) 


3 • / • 

ui = —, iro, 


Wi 


3 

U 3 = 


1 


in the domain of definition of ujj. 


(ii) In the rest of the section we will study the effect of the introduction of some 
possible regularizations, namely two of them 


• du* = CJ5 dv. 


• du* = LO 3 W4 015 du. 

Again, we have to keep in mind that with the notation used in the above regu¬ 
larizations, the new variable v* is different from one case to the other. 


7.1 The dv* /dv = (u 5 regularization. 

Then, associated to the ratios Ui , if we carry out the regularization 


du = u,5 du*. 


( 94 ) 


we have the following regularized differential system 


dui 
du* 
d u 2 

du* 

du 3 

dv* 

di/4 

dv* 


Cfu 2 U3U4, 
C'l U1U3U4, 

C 3 U\U-2U4. 

C\ uiu 2 u 3 , 


( 95 ) 


with IC Ui( 0 ) = i = 1 ,..., 4 . 

Thus, dividing the integral oiu;| — = C\ by we write: u\ = (04 — 

a^u\)/C\. Then, we obtain 


v = 



?i2 [ui(u*)] 2 du*, 


( 96 ) 


where n 2 = 0:5/01 and rti(u*) is a function solution of the system ( 95 ); quadra¬ 
ture which will be solved numerically. As we know that the solution of ( 95 ) can 
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be obtained by undetermined coefficients, making use of the A-Mahler functions 
defined by the system (76), but in the variable v*. In other words, the previous 
form of the solution represents an alternative to the use of hyperelliptic integrals 
for solving ([Tj) for N = 5. Or, in a more precise form, we have separated geom¬ 
etry from dynamics. The trajectory is expressed by Jacobi or Mahler functions, 
meanwhile the quadrature of the parametrization (96) will lead generically to a 
hyperellictic integral. 


7.2 The dv*/dv = U 3 U 4 UJ 5 regularization. 

Let us consider again the system 5-EES ( ]92| ). Now we try the regularization 

du* 

=w 3 w 4 u; 5 . (97) 

dv 

in a domain where U 3 UJ 4 U 5 7 ^ 0. This means that the system reduces to 


dcui 

du* 


— U>2, 


d C 02 
dv* 


= a 2 uj 1 , 


(98) 


and three quadratures associated to w 3 , 0 J 4 and 0 J 5 . In fact, they are not needed 
because the integrals allow to write ujf = C{ — aiiJ\ , (i = 3,4, 5). Remember that 
C\ are constants, functions of the initial conditions. 

Assuming the bounded case we can always choose, by scaling and transfor¬ 
mation of functions, that a 4 = 1, 012 = — 1- In other words we have 


oji(v*) = sinu*, 0J2{v*) = cosu*, 


(99) 


Then, the quadrature (97), taking into account the previous mentioned integrals, 
we have 

A„ = [ dt, ‘ (100) 




where /% and A are functions of C\ and 07 . This lead us, in the generic case, to 
a hyperelliptic quadrature. 


Dealing with the 5-Mahler System. In what follows we choose as the basic 
system in N = 5 a Mahler type system 

Uj[ = OJ2UJ3UJ4 UJ5, 

J 2 = -uji w 3 w 4 w 5 , 

W3 = — m cjio; 2 w 4 W5, (101) 

w 4 = -nwiw 2 w 3 W5, 

CJ5 = — pU\ U2UJ3U4, 

with initial conditions ( 0 , 1 , 1 , 1 , 1 ). 
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Moreover, apart from adjusting coefficients, an alternative form of dealing 
with (100) is to make a change of variable sinx* = x. Then, the corresponding 


new expression for the regularization is given by 


Adu = 


dx 


\J{\ — x 2 )(l — mx 2 )( 1 — nx 2 )(l — px 2 ) 


( 102 ) 


Denoting 


w = A v 


we define by Amg (generalized amplitude) the inverse function 

u* = Amg (w,p,m,n). (103) 

Then, by analogy with the notation introduced in lower dimensions, we propose 
to write 

sin v* = sin Amg (w;p,n,m) = Sng (w;p,n,m) (104) 

In other words, we define Sng 

x = x{w\p , n, m) = Sng(ie;p, n, rn) (105) 

as the three-parameter function solution of the differential equation 
/ Ht \ 2 

(-—) = (1 — x 2 )(l — px 2 )(l — nx 2 )(l — mx 2 ). (106) 

V duj J 

In the rest of this paper we restrict ourselves to the domain of parameters A = 
{(p,n,m) G [0,1] x [0,1] x [0,1]}. 

Then, associated with Sng we introduce the following functions 


Cng (w;p, n,m) = ±\J 1 — Sng 2 (w,p, n, m), 

Dng (w,p,n,m) = \J 1 — m Sng 2 (tc;p, n, m), 

Fng (w;p,n,m) = \J 1 — n Sng 2 (w, p, n , rn), 

Hng (w,p,n,m) = \J 1 — p Sng 2 (w,p, n, m). 

To simplify the notation, we will write in some expressions 

Sng(u;; p, n, m) = Sng, Cng(u;; p, n, m) = Cng, 
Dng(tc; p, n, m) = Dng, Fng(ie; p, n, m ) = Fng, 
Hng(u;; p, n, m) = Hng. 


(107) 
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Then, we write again (101) as the following IVP 


dSng 

du; 

dCng 

du; 

dDng 
dtc 
dFng 
d w 
dHng 
d w 


= Cng Dng Fng Hng, 

= —Sng Dng Fng Hng 
= —m Sng Cng Fng Hng 
= —n Sng Cng Dng Hng, 
= —p Sng Cng Dng Fng 


( 108 ) 


with initial conditions (0,1,1,1,1). 
grals take the following form 


Note that in agreement with (107), the inte- 


Cng 2 + Sng 2 = 1, Dng 2 + mSng 2 = 1, 
Fng 2 + n Sng 2 = 1, Hng 2 +pSng 2 = 1. 


(109) 


We are not going to deal with the generic study of our system (108). It is out 
of the scope of this paper. In the last Section we will restrict to analyze some 
particular cases 


8 N = 5 : Some particular cases 

Like in previous dimensions, we consider two particular cases 


8.1 The case p = 0. 


Now, according to (107), we have Hng 
studied case: 4-Mahler system. 


This corresponds to the previous 


8.2 The case p = n. 


As we have just pointed out, a particular case of (100) we will consider now two 
of the fa equals. According to the notation introduced, we write 

dd 


Xv = 


'o (1 — nsin 2 d) \ZY--rnsIrPd 


( 110 ) 


Remark 8.1 In relation to the quadrature ( fii0| ) the reader will remember that 
this is precisely the Legendre third elliptic integral n(7*; m, n). Thus, for the 

2 Dealing with the search of fast numerical algorithms for the computation of the third elliptic 
integral Fukushima HE] singles out in a recent paper that by a number of transformations the 
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particular cases n = 0 and n = m, we encounter the other Legendre elliptic 
integrals: 


F(ip, m) = 




dd 


\J\ — m sin 2 d 
E(tp,m) = / \J \ — m sin 2 d dd, 

Jo 

= (1 — m) II( 93 , m, m) + m 


= 11(93, 0 , m ) 


sin(2</3) 


2 ^ 


m sin 93 


Denoting 


w = Xv, 


we define as Amg (generalized amplitude) the inverse function 


v* = Amg (w,n,n,m). 


(Ill) 


Then, by analogy with the notation introduced in lower dimensions, we propose 
to write 


sin v* = sin Amg(tu; n, n, m) = Sng(ic; n, m ) (H2) 

For later use, we also include here the expression for our particular case of ( |102 ) 

dx 


d w = 


(1 — nx 2 )\J{l — x 2 )(l — mx 2 ) 


(113) 


From our initial conditions we have Hng = Fng. Then, from (101) we immediately 
obtain that (^4 = ^ 5 , and that these functions satisfy the following IVP 


^ f = CngDngFng 2 , 
aw 

^ y — = — SngDngFng 2 
aw 

dDn g o ri t? 2 

—-- = — m brig Cng h ng 

aw 

——— = —nSng CngDngFng, 


con las condiciones iniciales ( 0 , 1 , 1 , 1 ). 

Again by a regularization w —> v given by 


dx 
d w 


Fng 


(114) 


(115) 


domain of n and m are reduced as 

„ , ,— m 

0 < m < 1, — \Jm < n < - , 

1 +VT -m 

This fact has to be in mind (incluir grfico de este dominio) in order to study the Wj thinking on 
applications to those integrals... 
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transforms (114) in a regularized system which is a 4-Mahler system in the new 
variable. 

After we have solved the regularized system, we still need to compute the 
d w = I (116) 


quadrature associated to the differential relation (115). Explicitly we have 

dv 


yjl — nSng 2 (w(v)) 


We will give details of this process, both from the analytical and numerical point 
of view, in a forthcoming paper. 


9 On the application to the free rigid body 

We will apply what we have presented in previous sections to the description of 
the solution of the free rigid body. We will do that formulating the system in 
symplectic Andoyer variables. 


9.1 The solution in Andoyer variables 


Let us consider the Hamiltonian of the free rigid body expressed in Andoyer’s 
variables (A, /i, v, A, M, N) which takes the form 

T~L = ^(ai sin 2 v + 02 cos 2 v)(M 2 — N 2 ) + ^-IV 2 , (117) 


where ( 01 , 02 , 03 ) = (1/A,1/B,1/C) with ( A,B,C) the principal moments of 
inertia. Note that in applications we will study the influence of B, which will be 
taken as physical parameter A < B < C, join with C < A + B. The differential 
system is given by three equations 


d v 
dt 
d N 
dt 

dfi 

dt 


&H 

~dN 

&H 

dv 

d H 
dM 


= N{a 3 — 01 sin 2 v — 02 cos 2 v), 

= (02 — oi)(M 2 — N 2 ) sin v cos v, 
= M (01 sin 2 v + a ,2 cos 2 v), 


(118) 

(119) 

( 120 ) 


and the other three (A, A, M) which are integrals. Usually we integrate first the 
system defined by N and v. More precisely, we solve the Euler system, associated 
with those variables. Then, the functions solution N[t) and v(t) are given making 
use of the Jacobi elliptic functions 


sinzz(t) = 


cn(st; m) 


\/l + n*sn 2 (s t; m)' 

sn(st; m) 


cos v(t) = v 7 1 + n* - _, 

yl + n*sn 2 (s t; m) 

N(t) = Rdn(st] m), 


( 121 ) 
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where 


R 2 


= M 2 


C( 1 - <L4) 
C - A 


(B-A)(dC-l) 

(C-B){l-dA) 


n = 


C(B-A) 

" “ A{C-By 
(C-B)(l-dA) 
ABC 


with d = 2 h/M 2 . 


Remark 9.1 The reader will notice that sin is(t) and cosz^(t) are Mahler func¬ 
tions. Moreover from (78\ ) we know that N(t) is a ratio of Mahler functions 

Finally, we obtain n(t) by means of a quadrature: 

H = M J (ai sin 2 u{t) + <22 cos 2 v(t )) d t (122) 

which is finally expressed by means of a linear function of time and the Legendre 
third elliptic integral. Integral whose solution Jacobi gave making use of his 
elliptic and related functions. 


9.2 On alternative approaches 


We proceed here in a different form than the previous Section |9.1| Making use 
of the Hamiltonian function, we may separate variables in the system defined by 
(118)-(119). More precisely, we denote 


a\ — a 2 ai — 02 

ni = — -, mi = - 

d — 02 03 — 02 


(123) 


and n = (d — 02 X 03 — 02 ), where we assume 03 X 02 and d X 02 ; (the case of 
equality has to be treated separately). Then, the equation (118) may be written 
in the form 

Mn dt = dU —— (124) 

\J (1 — ni sin 2 u){ 1 — m\ sin 2 u) 


where n\ = n±(d, 02 ) and m\ = m\{d , 02 ), that is to say, we may study the system 
under the influence of the intermediate moment of inertia and the value of the 
Hamiltonian, keeping fixed the other parameters. Again, we have to distinguish 
circulation and libration patterns, but we do not need to go into details of that 
procedure here. 


In a more detailed form, we see that from (40) and (121) we may write 


sin u(t) = A\ eng (w, n, m), 
cos v{t) = A 2 sng(wj; n, m), 

N = A 3 dng(w;n,m)/fng(w,n,m), 


where Ai are quantities depending on the previous constants. 
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The quadrature (122) of the Andoyer angle variable /jl, now takes the form: 


(i = M J (ai sin 2 v + a ,2 cos 2 v) d t 
= M (oisng 2 t + 02 cng 2 t) d t 


= Ma 2 t +a\ J sng 2 fdt 

Finally, from what we have seen in Sect. [8] we find that 

sin [i = Sng(ir; n, rn ), cos (jl = Cng(u;; n, m) 


(125) 


In other words, depending on the use of sng, etc. or Sng, etc. we reach the third 
Legendre elliptic integral in two different forms. Comparisons of the pros and 
cons of their use, versus the classic approach based on sn, etc. Jacobi functions, 
is in progress. 


10 Appendices 


Appendix A: On the ratios of Jacobi 0* functions as solutions of 3-EES. 

From Lawden [9] (Clip. 1) we borrow the following 3-EES differential systems 
satisfied by the ratios of the Jacobi 0* functions 


d_ 

dv 

d_ 

dv 

d_ 

du 



(126) 

(127) 

(128) 


etc. We hnd convenient to introduce the notation = 0j/0 ? ; and the reparametriza- 
tion v —> t given by dr = ^2K/tt dv, with x'- = dxij/dr. Thus, taking into 
account the values of 0j(O), where k 2 = m, k 2 + k ' 2 = 1 and K(m) is the the 
complete Legendre first elliptic integral, we write those IVP systems as follows. 
Note that, as was pointed out in Crespo and Ferrer [3j, considering the sign of 
the coefficients, we may distinguish 


• Two bounded systems: 


3^41 — k X42 X43, 

X42 = — X41 X43, 

X43 = —k X 41 X 42 , ( 0 , \/k/k ', 1 /VW) 


and 


3-31 = ^ 32 ^ 34 , 

X32 = — k X31 X34,, 

X34 = k £31 x 3 2 , (0, Vk, Vkf) 
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Two unbounded systems: 


x' 2 i = kx 23 ^ 24 , 

3-23 = ^ X 21 %24i j 

®24 = X 21 ^23 j (0, l/\/fc, \fkfjk) 


kx 13x14, 

®12 ® 14 , , 

k'x 12 xi 3 , (1, V(&' + 1 )/k, V(k' + l)/k). 

Then, we may express those ratios as functions the Jacobi elliptic functions and 
their Glashier ratios. 

Appendix B: Transformations and addition formulas for Jacobi elliptic 
functions. 

For the benefit of the reader we bring here some well known transformations 
involving the elliptic modulus. They may be found in any handbook of elliptic 
functions (remember that, depending on the authors, two notations are used: 
‘modulus’ or ‘parameter’ related by k 2 = m, and their complementaries). Those 
formulas should be used for the reduction to the normal case of some of the 
particular cases mentioned along the paper. 

• Negative parameter 
Let m be a positive number and write 

m 1 u 

-’ ki = , V= N=- 

1 + m 1 + m yjpi 



and 


x'12 

x 13 

x 14 


Then, 


sn(it; —m) 
cn (u ; —m) 
dn(u; — m) 


\[W\ 


sn(u; p) 


dn(u; p ) : 
cn(u; p) 
dn(u; p) ’ 


1 


dn(v ; p) ’ 


Thus elliptic functions with negative parameter may be expressed by elliptic 
functions with a positive parameter. Note that 0 < p < 1. 

A final comment related to the complete elliptic integral of first kind is due 
here. Unlike Maple, the software Mathematica yields the following result 

d(f> ^ K f m 

\/\ — m sin 2 cj) \/l — m \ 171 ~ 1 
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for Vm < 1, instead of the expected result K(m). By applying the previous 
change (129), we have that, being m a positive number, 


K(-m) = —=L== K f yj—) = \/Mi k (m) (131) 

y/l + m \l + mj 

which is exactly the same result given by Mathematica for m < 0. 

• Reciprocal parameter 
Denoting now v = \fmu, we have 

sn (u ; m) = —==. sn(n ; m ~ 1 ), 

\Jm 

cn (u ; m) = dn(u ; m -1 ), 
dn(u ; m) = cn(u ; m _1 ). 

This is Jacobi’s real transformation. If m > 1, then m~ l < 1, thus elliptic func¬ 
tions whose parameter is greater than 1 are related to the ones whose parameter 
is less than 1. In short there is no loss of generality assuming 0 < m < 1. 

• Decrease of parameter 


d = 


/ I ~ Vrnj \ 2 

VI + yfm{) 


u 

1 + \fd 


(132) 


sn(u; m) 
cn (it; m) 


dn(u; m) 


(1 + \//l)sn(u ; p) 

1 + \/~R sn 2 (v ; jLt) ’ 

cn(n ; p) dn(n ; p) 

1 + \[d sn 2 (v ; p) ’ 

1 - y/~psn 2 {v ; p) 

1 + ^/p sn 2 (u; p)' 


This is Gauss transformation or the descending Landen transformation, which 
makes elliptic functions to depend on functions with a smaller parameter. 

Note that, making use of the double angle, we may also write 


dn(u; m) 


^fp cn(2w ; p) + dn(2u ; p) 

1 + yfd 


(133) 


There are analogous expressions for the increase of parameter. For a recent study 
where generalized for mules are given, see [8J. 

• Addition formulae 

Complementing previous transformations, we collect also here the addition for- 
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mulae 


sn(a + /3) 
cn (a + f3) 
dn(a + j3) 


sn a cn /3 dn [3 + sn ft cn a dn a 
1 — m sn 2 a sn 2 ft 
cn a cn ft — sn a sn ft dna dn /3 
1 — m sn 2 a sn 2 /J 
dn a dn ft - msnasn^cna cn /3 
1 — m sn 2 a sn 2 /3 


which we have generalized for the new functions; more precisely this has been 
done for the 4-EES Mahler system. 
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